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SUMMARY 


Multigrid methods are good candidates for the resolution of the system arising in Numerical 
Fluid Dynamics. However, the question is to know if those algorithms which are efficient for the 
Poisson equation on structured meshes will still apply well to the Euler and Navier-Stokes equations 
on unstructured meshes. The study of elliptic problems leads us to define the conditions where a Full 
Multigrid strategy has 0(N ) complexity. The aim of this paper is to build a comparison between the 
elliptic theory and practical CFD problems. 

First, as an introduction, we will recall some basic definitions and theorems applied to a model 
problem. The goal of this section is to point out the different properties that we need to produce 
an FMG algorithm with O(N) complexity. Then, we will show how we can apply this theory to the 
fluid dynamics equations such as Euler and Navier-Stokes equations. At last, we present some results 
which are 2nd-order accurate and some explanations about the behaviour of the FMG process. 


INTRODUCTION 


One first important element is the mesh independent convergence speed. Hackbush, in [1] for 
example, proposes a demonstration of this property. It is done in the special case of an elliptic 
problem on structured nested meshes. We want to evaluate the properties that we must keep in order 
to get the mesh independent convergence speed when we use unstructured non-embedded meshes. 

The problem to be solved is the following: 

{ Att = / onfl convex polygonal domain 

0 ^ 

u E H?(Q) and f E L 2 (D) 

The discretization is a usual linear Pl-Galerkin finite element. Thus, we get a discrete space 
whose dimension is equal to N h (number of nodes), and where the subscript h indicates the mesh 
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size. The resulting problem consists now in solving the linear system: 


Ah Uh — fh (2) 

We may evaluate the discretization error thanks to the Aubin-Nitsche’s theorem and usual regularity: 

\\u — Uh\\L 2 < C 2 h 2 \\u \\ H 2 < Cz h 2 H/Hn 2 (3) 

The linear system (2) may incur a lot of CPU time because of its size (large number of nodes). The 
idea is then to use a second finite element subspace Tin whose dimension N H is less than the previous 
one (usually H = 2 h). We have then the following relationship between both spaces (also called 
grids): 

At 1 


X H 
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1R Nh - 

-> R n " 
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-*• ffi Nh 


A~ l 
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where P and R are linear interpolations (transfer operators). The iterative process can be written 
as: 


< +1 = M h < + N h h where: M h = S% (I - P A~ H l R A h ) , S h = I-uD~ l A h 

(S defines the basic iterative smoother, v\ and 1/2 the number of pre- and post-relaxations). Such 

a process converges if ||M/,|| < 1. A very important property of this kind of method is that the 
convergence is independent of the mesh size. In order to simplify the notations (and the study) we 
rewrite M h as the following ideal- 2-grid operator [2]: 

Mh = (V - P An K) W 

The norms of both factors of the right hand side of the equation (4) will determine the norm of Mh : 

• The smoothing property: 

MfcSfcll < 1/h 2 tj(v) 

lim ri{v) = 0 ' ' 

depends a lot on the basic smoothing process, and, we will not give any details. 

• The approximation property is: 

IK 1 -PA„ l R || = 0(h?) (6) 

Let us focus on (6): it takes into account the transfer operators, and overall, represents the 
difference that exists between the solution on the fine grid and the solution on the coarse grid. 

An MG scheme that exhibits these properties will result in a convergence speed that is independent 
of the mesh size: 

Vp 3 i/(p) such that || M h v h - u fc || < p\\v h - u h || 

We may notice that demonstrating the approximation property leads to the evaluation of the following 
quantity: 

|| phPAjfRrh - A -1 1| 
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For nested meshes, thanks to (3), one can easily derive this from the following equality: 


PhP = Ph 

On the other hand, for unnested meshes, p^P is not equal to pn and this evaluation is more difficult. 
Actually, it is the same as evaluating the difference between two interpolated solutions. Zhang [3], 
thanks to Bank-Dupont’s theory, proposes such an evaluation. 

Remark: The multigrid iterative V-cycle algorithm can easily be deduced from the previous 2-grid 
algorithm recursively. It maintains the convergence speed independently of the mesh size and has an 
0(N log N) complexity. 

In order to apply the previous result of convergence, we propose to use a Full Multigrid (FMG) 
strategy (proposed by Brandt in [4]). A well known result is the one given by Hackbush in [1], which 
is given below: 

Theorem: We note: k the consistency order, Ci = max (hk-i/hk) K the ratio of accuracy 

1 < k < l 

between two solutions, S the reduction factor of the MG process, * the number of MG iterations 
applied to reach the solution u \ . If Uk is the solution of the discrete problem, we have: 

11*4-1**11 < ciCtK with ci = l _ S c jS , 

Assuming that C% = 2*, we deduce the exact number of cycles in each FMG phase to solve the 
lst-order problem ( S' < 1/4) and the 2nd-order problem (S'* < 1/8). The number of cycles i in 
each phase is constant, which leads to an algorithm that has O(N) complexity. 

Once again, the relative interpolation error [1] conditions the quality of the initialization in each 
phase. Thus, in order to stay close to the ideal scheme (where the different subspaces are nested), we 
propose to build meshes where: 

• The mesh size ratio is close to 2, 

• The triangles aspect ratio is locally comparable in the whole domain. 

We have thus identified the different necessary ingredients to build an algorithm having O(N) com- 
plexity: 


1. A sequence of grids, 

2. A basic smoother (ex: Jacobi, Gauss-Seidel), 

3. Intergrid transfer operators (ex: linear interpolations), 

4. MG algorithm (ex: V-cycle, W-cycle), 

5. FMG strategy. 

We may now apply it to more complex fluid dynamics problems such as the resolution of the 
Navier-Stokes/Euler equations. 
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FLUID DYNAMICS APPLICATIONS 


>C 


We recall first the formulation of the steady Navier-Stokes equations: 
dF(W) dG(W ) 1 (dR{W) dS(W)\ 


F(W) = 


R{W) = 


/ pu \ 


( pv > 

pu 2 4 -p 

G{W) = 

puv 

2 1 

puv 


pv* + p 

( E+p)u J 


( E+p)v 
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•xy 


'ynde 

UT XX + vr xy 


S{W) = 


= ( $$ ) 
0 \ 


I xy 

T yy 

ur xy + vr yy +^-— J 


( 7 ) 


V = (7 - !) ( E ~ \ P (« 2 + ^ 2 )) . e = C V T = j~ \{u 2 + v 2 ) , 

where 7 = 1.4 is the ratio of specific heats, T is the temperature, \i and k are the normalized viscosity 
and thermal < 
are given by: 


and thermal conductivity coefficients. The components of the Cauchy stress tensor r XI , r xy and r yy 


2 ( n du dv\ 2 / dv du\ ( du dv 

Txx ~ 3^ \ ~dy) ’ Tyy = 3 M V^ %~di) ' Txv = fi [d^ + di 

Re = poUoL 0 /fj, 0 is the Reynolds number and Pr — hqC p /kq is the Prandtl number, where p 0 , Uo, L 0 
and Hq denote respectively the characteristic density, velocity, length and diffusivity of the considered 
flow. It is easily seen that, if the right-hand side is equal to zero, then we recover the Euler equations. 
The inlet conditions are defined by the farfield flow. For Euler flows, we impose the slip condition 
on the wall (V.n = 0), and for Navier-Stokes flows, on the wall, the no-slip condition (V = 0) and 
the isothermal condition (T = TJ,). The discretization is given by a mixed FEM/FVM formulation 
[5], where the mesh is a finite-element type (triangles), on which we construct control-cells (FVM) in 
order to solve the variational formulation of the equations, such as, for the Euler flows: 



E 


f i' i .lF(W)d 

JdQj 


+ I 

JTndC, 


v i .F(W)da = 
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j dxdy 


( 8 ) 


The computation of the fluxes, appearing in (8), between two cells, is managed by Roe’s numerical 
flux vector splitting in the domain, and by the Steger- Warming numerical flux vector splitting for the 
farfield boundaries. The 2nd-order accurate scheme is obtained by the use of the MUSCL method 
developed by van Leer [6]. We solve the discrete equations with non-linear relaxation algorithms [6], 
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namely here the multistage Jacobi algorithm [2, 7]: 

wj 0) = Wf 

For ks = 1 to nstage 

( 9 ) 

wf s+1) = wj 0) -to c ke J2 Wj » w j k9) > ^.j) 

ieK(j) 

p^Oc+1 __ yyi n3 t a 9 e ) 

Let us now look at the meshes. We start with an initial given fine mesh Fig.l.a. Finer meshes are 


Pki = 



a. Initial 800 nodes 



b. 3114 nodes 


c. Finest 12284 nodes 


Figure 1: NACA0012 fine meshes. 


obtained by triangle subdivision (Fig.l.b,.c). Then, we use a coarsening algorithm due to Guillard 
[8] to build coarser meshes, from the initial one. This produces a sequence of node-embedded meshes 
(Fig.2.a,.b,.c). We get 6 meshes for the NACA0012 profile, where the finest has 12284 nodes and the 
coarsest 19 nodes. This method allows us to keep the mesh size ratio close to 2, and a comparable 
local mesh aspect ratio. The intergrid transfer operators [9] are linear interpolations, concerning the 
variables and the corrections, and linear distributions, concerning the residuals. The MG algorithm 
will be the W-cycle, because it is the natural extension of the ideal-2-grid scheme. Furthermore, there 
exist several ways to obtain 2nd-order accurate solutions: 

• Mavriplis [10] uses an FMG algorithm, where lst-order accurate solutions are computed on the 
coarse levels, and 2nd-order accurate on the finest. Some experiments with our upwind schemes 
showed us that the convergence speed is hardly independent of the mesh size. 

• Hemker-Koren [11] propose to get a lst-order solution with an FMG strategy and then to 
compute a certain number of DeCV-cycles: They use the Defect Correction (DeC) algorithm 
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a. 223 nodes 




Figure 2: NACA0012 coarse meshes. 


[12], in order to solve the 2nd-order accurate following problem: TiiW) = S. It is written: 


FiiW 1 ) = S , 

jF 1 (W Ar+1 ) = JF 1 (lF JV )-^ 2 (VF JV ) + 5, N = 1,2,... 


( 10 ) 


Actually, we define the DeCV-cycling method where the lst-order problem in (10) is approxi- 
mately solved with one V-cycle. However, we do not know how many DeCV-cycles are to be 
performed and we lose the O(N) complexity. 

We propose here two different methods in order to obtain 2nd-order accurate solutions with an 0(N ) 
complexity algorithm [1]. 

• FMDeCV is an FMG strategy where we use DeCV-cycles in each phase, with two Jacobi sweeps 
per level (FMDeCV-2RKl). 

• FMG2 is an FMG strategy where we use on each level of the different phases the good damping 
properties of the multistage schemes (see [13]) for smoothing directly the second order accurate 
problem. 

A result of convergence of the DeCV method is given in [14] and assures that DeCV-cycling has a 
convergence speed independent of the mesh size: 

\\DeCVUfr — < S 2 U^|| , S 2 = Si 4- SiSoeC + $DeC < 1 (11) 

where DeCVu £ is the a — th iterate of the DeCV-cycling and u\ is the solution of the 2nd-order 
accurate problem. 


Remark: Desideri-Hemker in [12] show that the convergence speed SoeC of the DeC process is at 
least equal to 1/2. From (11), we should use MG lst-order accurate algorithms whose convergence 
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IILog (RES)II 


speed Si (independent of the mesh size) is less than 1/3. Actually, the experiments showed us that 
S\ > 1/3 does not induce difficulty. 


SECOND ORDER ACCURATE RESULTS 
Euler Flows around a NACA0012 profile 







Figure 3: Euler FMG2 phases, isomach contours, M x = 0.9, a — 0°. 

The test-case depicted in Fig. 3 to Fig. 5 is defined by a farfield Mach number equal to 0.9, and 
a zero angle of attack. The smoother is the (4 stage) RKJ one, whose coefficients (oi 0.14, 0:2 
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Contours interval: 0.50000E-01 

Min./ Max. Mach 0.12322, 1.4793 

b. 3114 nodes c. 12284 nodes 

Figure 4: Euler FMG2 phases, isomach contours, = 0.9, a = 0°. 

0.2939,03 = 0.5252,a 4 = 1) are due to van Leer [15], and defines the FMG2 strategy. In Fig.3.a, 
we present the convergence histories of the logarithm of the residual versus the number of cycles in 
each phase. The convergence is estimated as obtained when the residual decrease reaches 10~ 6 . The 
first history is a 1-grid convergence on the coarsest mesh (19 nodes), the last (4-grid scheme) on the 
800 node mesh. At the end of the convergence of each phase we produce an initialization of the 
next phase by interpolating the solution on the next finer mesh. We may notice that the different 
convergence histories tend to be straight lines with the same value of slope: this allows us to say that 
the convergence speed is independent of the mesh-size and that we may use an FMG strategy. In 
Fig.3.b the residual convergence histories of the last 5 phases are depicted (the first one is a 1-grid 
convergence history with a residual decrease equal to 10~ 6 ). In order to neglect oscillatory non-linear 
phenomena we choose to impose a residual decrease equal to 1/8 (and not a defined number of cycles): 
the FMG convergence histories follow exactly the peaks of the (corresponding) phases on Fig.3.a, and, 
the solutions are not either changed. A solution on the finest grid (Fig.l.c) is reached after only 5 
cycles (77 WU, 673 s on Convex C210 with non vectorized software). In Fig.3 and Fig.4 we show the 
different solutions obtained at the end of each phase: they are non-symmetrical solutions (Fig.3), due 
to the fact that the coarse meshes are non-symmetrical. However, this phenomenon vanishes when 
the different meshes become symmetrical (Fig.4). Another important remark is that the solution 
between the finest mesh and the next coarser one does not vary much: we may say that we get a 
nearly converged solution on the 3114 node mesh. In order to verify the previous assumption, we 
compare the FMG solution of Fig.5 with the solution obtained after a 10 -6 residual decrease: the 
Mach number extrema are approximately the same although the isomach fines are a little bit more 
oscillatory for the FMG solution. 



Contours interval: 0.50000E-01 

Min./ Max. Mach 0.45E-01, 1.519^ 



Contours interval: 0.50000E-01 

Min./ Max. Mach 0.22627,1.4630 


a. 800 nodes 
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Contours interval: 0.50000E-01 

Min./ Max. Mach 0.45E-01, 1.5196 


a. FMG2 



b. Converged (RES = 10 6 ) 


Figure 5: Euler FMG2 solutions, isomach contours, M ro = 0.9, a — 0°. 


Navier-Stokes Flows around a NACA0012 profile 


The first test-case is defined by a farfield Mach number equal to 0.8, an angle of attack equal to 10 
degrees and a Reynolds number equal to 73. We use here an FMDeCV-2RKl strategy (justified by 
the mesh independent convergence of Fig. 6. a). The solution (Fig.6.c) is obtained on the finest mesh 
after 4 cycles that represent 42 WU computation and a total CPU time equal to 537 s. In Fig.6.d we 
illustrate the behavior of the pressure lift, CL (drag, CD) coefficient with a solid (dash) fine, versus 
the number of finest-grid iterations, up to a residual decrease on the finest grid of 10 12 . The points 
(cross) represent these coefficients during the FMG process, thus up to a 0.125 residual decrease. We 
can notice that the value of each of them is almost obtained at the end of the FMG phase (the error 
is equal to 64. 10~ 5 for the CL coefficient and to 16. 10 -5 for the CD coefficient). 

The second test-case, presented in Fig.7, is defined by a farfield Mach number equal to 2, an angle 
of attack equal to 10 degrees, and a Reynolds number equal to 106. This time we use an FMG2 
strategy (more robust than FMDeCV-2RKl that needs TVD limitation and implies that the residual 
stalls from the value of 10 -3 ); this produces a solution after 7 cycles (Fig.7.c, 109 WU, 1862 s), and 
we can make the same remarks as in the previous test-cases. The next coarser mesh results in CL 
and CD values within 1% of their final values which are obtained after 7 cycles on the finest mesh 
with a related error respectively of 2.10 -6 and 2.10 -5 (Fig.7.d), 

The last test-case shows us one limitation of this method. It is defined by a farfield Mach number 
equal to 0.8, an angle of attack equal to 10 degrees and a Reynolds number equal to 500. Here again, 
we use an FMDeCV-2RKl strategy (Fig.8). We may note, once again, that the FMG convergence 
histories (Fig.8.b) look like the corresponding peaks on Fig.8.a, thus we think that the solution does 
not vary between a 10' 1 residual decrease and a 10 -6 one. However, on Fig.9.a we note that the 
isomach lines are deformed up until the last solution (Fig.9.c) which presents two bumps. Actually, 
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Figure 7: Navier-Stokes FMG2 solution, M x — 2, a — 10°, He — 106. 
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LAST 5 PHASES RES(x) = 0. 125 



Contours interval: 0.25000E-01 Contours interval: 0.25000E-01 Contours interval: 0.25000E-01 


Min./ Max. Mach 0. ,0.79407 Min./ Max. Mach 0. , 1.0496 Min./ Max. Mach 0. , 1.0901 


c. 19 nodes 


d. 67 nodes 


e. 223 nodes 


Figure 8: Navier-Stokes FMDeCV-2RKl phases, isomach contours, M x , = 0.8, a = 10°, Re = 500. 
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since the solution differs a lot between the one obtained on the finest mesh and the other on the 
next coarser mesh (Fig.9.c and .b), we may say that we did not obtain a grid-converged solution. 
Furthermore, these two humps may be justified because the meshes that we use (especially the coarse 
ones) are not adapted for such a viscous flow and may not capture the boundary layer. Our assumption 
is confirmed by the CL and CD behaviors (Fig.lO.b), where we note an error for the CL coefficient 
equal to 5.1CT 2 and for the CD coefficient equal to 8.10 -3 . Moreover, we get a solution (Fig. 10. a) 
after 105 cycles (1107 WU, 14112 s), for a residual decrease of 10 -12 and which tends to confirm that 
the FMG solution is not a steady solution. 



Figure 11: Convergences 
CONCLUSION 


We want to point out that the use of FMG2 or FMDeCV strategy allowed us to get 2nd-order 
accurate solutions in most of cases with a limited number of operations (O(JV) complexity) . FMDeCV- 
2RK1 is more adapted to smooth problems and costs half as much as FMG2. Furthermore, we had 
to use an entropy correction technique [16] to get the above results, and occasionally a 10~ 12 residual 
decrease required to increase the value of this correction to prevent the residual from stalling or to 
improve robustness on the finest grid (this increased slightly the number of cycles in each phase 
without changing the solution greatly). As depicted in Fig. 11, these types of computations do not 
induce any stall during the convergence. The main two difficulties, non-embedded meshes and the 
requirement of 2nd-order accuracy, were remedied, respectively, by using a coarsening algorithm based 
on a Voronoi technique, and basic iteration techniques that were sufficient smoothers. The difficulty 
encountered in using an FMG strategy, with our meshes, increased as the Reynolds number was raised. 
Actually, it is obvious that our meshes are not adapted to these computations and that boundary 
layer problems will need the production of stretched meshes and different specialized smoothers, as 
suggested in [14]. 
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